Lotka-Volterra equations, or predator-prey equations, are a pair of first-order nonlinear differential equations describing the interaction between a food source and its consumers. If the prey is represented by the variable x and the predators by y, the equations are

$\begin{array}{l}\frac{dx}{dt}=ax-bxy\\ \frac{dy}{dt}=-cy+dxy\end{array}$

The first term on each right-hand side represents population growth in the absence of interaction: exponential increase for the prey and exponential decrease for the predators. The nonlinear terms describe the result of predation, with the prey decreasing in number and the predators increasing.

The equations are written appear to have four independent parameters, but with the scaled variables

$x=\frac{c}{d}Xy=\frac{c}{b}Yt=\frac{1}{c}T$

can be written as a system with a single independent parameter:

$\begin{array}{l}\frac{dX}{dT}=\alpha X-XY\\ \frac{dY}{dT}=-Y+XY\end{array}\alpha =\frac{a}{c}$

It is now straightforward to integrate these equations for a variety of initial conditions and parameter. The prey is in green and the predators in red for the obvious connotations:

The exponential increase of the prey can be seen by setting the initial predator population to zero, and the exponential decrease of predators by setting the initial prey population to zero.

The critical points of the nonlinear system are found by setting the right-hand sides of the differential equations equal to zero and solving for variable values. The two critical points for the scaled pair of equations are

$[0,0]\mathrm{and}[1,\alpha ]$

Linearizing the pair of differential equation near the first critical point

$[X,Y]\sim [{\epsilon }_X,{\epsilon }_Y]\to \frac{d}{dt}\left[\begin{array}{c}{\epsilon }_X\\ {\epsilon }_Y\end{array}\right]=\left[\begin{array}{cc}\alpha & 0\\ 0& -1\end{array}\right]\left[\begin{array}{c}{\epsilon }_X\\ {\epsilon }_Y\end{array}\right]$

produces a coefficient matrix with two real eigenvalues of different sign, so that this critical point is a saddle point. Linearizing the pair of differential equation near the second critical point

$[X,Y]\sim [1+{\epsilon }_X,\alpha +{\epsilon }_Y]\to \frac{d}{dt}\left[\begin{array}{c}{\epsilon }_X\\ {\epsilon }_Y\end{array}\right]=\left[\begin{array}{cc}0& -1\\ \alpha & 0\end{array}\right]\left[\begin{array}{c}{\epsilon }_X\\ {\epsilon }_Y\end{array}\right]$

produces a coefficient matrix with pure imaginary eigenvalues, so that this critical point is a center. Phase-space orbits move around this center in a counterclockwise direction, since ${\stackrel{·}{\epsilon }}_Y>0$ for ${\epsilon }_X>0$ .

To display orbits in phase space, an energy integral is needed and easily found for this relatively simple system. Eliminating the differential with respect to time in the pair of equations gives

$\begin{array}{c}\frac{dX}{X(\alpha -Y)}=\frac{dY}{Y(X-1)}\\ \frac{X-1}{X}dX=\frac{\alpha -Y}{Y}dY\\ \\ X-lnX+Y-\alpha lnY=E\end{array}$

This energy integral can be rearranged to express one variable as a function of the other in terms of the Lambert W function. This function satisfies the implicit equation

$W\left(z\right)e^{W\left(z\right)}=z$

To determine $X\left(Y\right)$ multiply both sides of the energy integral by minus one, exponentiate and multiply again by the same factor:

$\begin{array}{c}-[X-lnX=E-Y+\alpha lnY]\\ -[X{e}^{-X}=Y^{-\alpha }{e}^{Y-E}]\\ X=-W[-Y^{-\alpha }{e}^{Y-E}]\end{array}$

For $Y\left(X\right)$ the multiplicative factor includes the inverse of the parameter:

$\begin{array}{c}-\frac{1}{\alpha }[Y-\alpha lnY=E-X+lnX]\\ -\frac{1}{\alpha }[Y{e}^{-Y/\alpha }=X^{-1/\alpha }{e}^{(X-E)/\alpha }]\\ Y=-\alpha W[-\frac{1}{\alpha }{X}^{-1/\alpha }{e}^{(X-E)/\alpha }]\end{array}$

Since the arguments of both functions are negative, the expressions are located in the double-valued region of the Lambert function, which looks like this:

One part of each phase-space orbit will come from the principal branch of the function, with the other part from the other real branch. Since this second branch has larger absolute values, the upper part of each orbit in phase space will come from it.

To determine the limits of variables in phase space, differentiate the energy integral implicitly with respect to each variable:

$\begin{array}{c}1-\frac{1}{X}+(1-\frac{\alpha }{Y})\frac{dY}{dX}=0\\ \frac{dY}{dX}=0\to X=1\end{array}\begin{array}{c}(1-\frac{1}{X})\frac{dX}{dY}+1-\frac{\alpha }{Y}+=0\\ \frac{dX}{dY}=0\to Y=\alpha \end{array}$

The limits of an orbit for a given variable occur when the other has the value of the center critical point. For $Y=\alpha$ these limits are

$X_{\min }=-W_0[-\frac{e^{\alpha -E}}{{\alpha }^{\alpha }}]X_{\max }=-W_{-1}[-\frac{e^{\alpha -E}}{{\alpha }^{\alpha }}]$

and for $X=1$ the limits are

$Y_{\min }=-\alpha W_0[-\frac{e^{(1-E)/\alpha }}{\alpha }]Y_{\max }=-\alpha W_{-1}[-\frac{e^{(1-E)/\alpha }}{\alpha }]$

where the two values come from the two real branches of the Lambert function as previously stated. The critical point is also the absolute minimum of the energy integral as a surface, as can be seen by forming its gradient

$\nabla E=[1-\frac{1}{X},1-\frac{\alpha }{Y}]$

which is identically zero at the critical point. A visualization of the energy surface, with the critical point in red, confirms this behavior:

With limits available for both variables as well as the energy, one can now visualize an orbit for arbitrary energy, along with one for $E=5$ in purple and the center critical point again in red:

These orbits are patently horizontal slices of the energy surface above.

There are several ways to express the period of orbits as summarized nicely in this review article, the most direct making use of expressions already determined. The differential of the scaled temporal variable can be equally expressed as

$dT=\frac{dX}{X(\alpha -Y)}=\frac{dY}{Y(X-1)}$

Integrating either expression along an entire orbit provides an integral representation of the period. First denote integrals on each constituent branch of the orbit as

$\begin{array}{l}{T}_{X_{\mathrm{lower}}}=\underset{X_{\min }}{\overset{X_{\max }}{\int }}\frac{dX}{\alpha X[1+W_0(-\frac{1}{\alpha }{X}^{-1/\alpha }{e}^{(X-E)/\alpha }\left)\right]}\\ T_{X_{\mathrm{upper}}}=-\underset{X_{\min }}{\overset{X_{\max }}{\int }}\frac{dX}{\alpha X[1+W_{-1}(-\frac{1}{\alpha }{X}^{-1/\alpha }{e}^{(X-E)/\alpha }\left)\right]}\\ T_{Y_{\mathrm{lower}}}=\underset{Y_{\min }}{\overset{Y_{\max }}{\int }}\frac{dY}{Y[1+W_0(-Y^{-\alpha }{e}^{Y-E}\left)\right]}\\ T_{Y_{\mathrm{upper}}}=-\underset{Y_{\min }}{\overset{Y_{\max }}{\int }}\frac{dY}{Y[1+W_{-1}(-Y^{-\alpha }{e}^{Y-E}\left)\right]}\end{array}$

where overall signs are chosen to produce positive values for counterclockwise motion. The period of the complete orbit is then

$T_{\mathrm{orbit}}=T_{X_{\mathrm{lower}}}+T_{X_{\mathrm{upper}}}=T_{Y_{\mathrm{lower}}}+T_{Y_{\mathrm{upper}}}$

These two expressions must give identical results by definition, so either can be used in practice. Taking the second as a bit more concise, the period as a function of energy is

The orbit period appears to be quite linear for small energies, but does not remain so for higher energies or low values of the single system parameter. Plotting the two parts of the period separately shows that they are not equal and their ratio is not constant, indicating that solutions to the pair of differential equations are not simply circular or Jacobi elliptic functions.

With an expression for period, numerical integration of the differential equations can be performed over just a single one. This facilitates visualization of the result in phase space with animation of the motion along the orbit:

Having integral representations for the two parts of the period allows one to write explicit expressions for the time to maxima and minima in each of the variables. Inspection of the phase-space orbits makes clear that the expression for a given initial variable is controlled by the value of the other initial variable relative to the center critical point, since this determines the branch of the relevant Lambert function. The time to maximum population of prey is thus

$T_{X_{\max }}=\{\begin{array}{cc}\underset{X_0}{\overset{X_{\max }}{\int }}\frac{dX}{\alpha X[1+W_0(-\frac{1}{\alpha }{X}^{-1/\alpha }{e}^{(X-E)/\alpha }\left)\right]}& ,Y_0\le \alpha \\ T_{X_{\mathrm{lower}}}+\underset{X_0}{\overset{X_{\min }}{\int }}\frac{dX}{\alpha X[1+W_{-1}(-\frac{1}{\alpha }{X}^{-1/\alpha }{e}^{(X-E)/\alpha }\left)\right]}& ,Y_0>\alpha \end{array}$

while the time to minimum prey is

$T_{X_{\min }}=\{\begin{array}{cc}{T}_{X_{\mathrm{upper}}}+\underset{X_0}{\overset{X_{\max }}{\int }}\frac{dX}{\alpha X[1+W_0(-\frac{1}{\alpha }{X}^{-1/\alpha }{e}^{(X-E)/\alpha }\left)\right]}& ,Y_0\le \alpha \\ \underset{X_0}{\overset{X_{\min }}{\int }}\frac{dX}{\alpha X[1+W_{-1}(-\frac{1}{\alpha }{X}^{-1/\alpha }{e}^{(X-E)/\alpha }\left)\right]}& ,Y_0>\alpha \end{array}$

Similarly, the time to maximum population of predators is

$T_{Y_{\max }}=\{\begin{array}{cc}{T}_{Y_{\mathrm{upper}}}-\underset{Y_0}{\overset{Y_{\min }}{\int }}\frac{dY}{Y[1+W_0(-Y^{-\alpha }{e}^{Y-E}\left)\right]}& ,X_0\le 1\\ -\underset{Y_0}{\overset{Y_{\max }}{\int }}\frac{dY}{Y[1+W_{-1}(-Y^{-\alpha }{e}^{Y-E}\left)\right]}& ,X_0>1\end{array}$

while the time to minimum predators is

$T_{Y_{\min }}=\{\begin{array}{cc}-\underset{Y_0}{\overset{Y_{\min }}{\int }}\frac{dY}{Y[1+W_0(-Y^{-\alpha }{e}^{Y-E}\left)\right]}& ,X_0\le 1\\ T_{Y_{\mathrm{lower}}}-\underset{Y_0}{\overset{Y_{\max }}{\int }}\frac{dY}{Y[1+W_{-1}(-Y^{-\alpha }{e}^{Y-E}\left)\right]}& ,X_0>1\end{array}$

These expressions can be verified visually by repeating the initial integration of the pair of differential equations for a single period and adding vertical lines at the indicated times:

Other quantities of interest are the times between maxima and minima in each of the populations. It might seem at first that expressions for such quantities will require four different options as combinations of the expressions just given, until one realizes that they remain in a single phase-space quadrant. This makes them easier to evaluate than anticipated:

$\begin{array}{l}{T}_{Y_{\max }}-T_{X_{\max }}=-\underset{1}{\overset{X_{\max }}{\int }}\frac{dX}{\alpha X[1+W_{-1}(-\frac{1}{\alpha }{X}^{-1/\alpha }{e}^{(X-E)/\alpha }\left)\right]}\\ \phantom{T_{Y_{\max }}-T_{X_{\max }}}=-\underset{\alpha }{\overset{Y_{\max }}{\int }}\frac{dY}{Y[1+W_{-1}(-Y^{-\alpha }{e}^{Y-E}\left)\right]}\end{array}$

$\begin{array}{l}{T}_{Y_{\min }}-T_{X_{\min }}=-\underset{1}{\overset{X_{\min }}{\int }}\frac{dX}{\alpha X[1+W_0(-\frac{1}{\alpha }{X}^{-1/\alpha }{e}^{(X-E)/\alpha }\left)\right]}\\ \phantom{T_{Y_{\min }}-T_{X_{\min }}}=-\underset{\alpha }{\overset{Y_{\min }}{\int }}\frac{dY}{Y[1+W_0(-Y^{-\alpha }{e}^{Y-E}\left)\right]}\end{array}$

Each of the pairs of integrals are easily seen to be equal: the integrands are equivalent by construction and the limits of integration follow from the corresponding change of variables. Again taking the simpler of the two representations, the temporal differences as functions of energy between maxima (in purple) and minima (in bright green) look like this:

The maxima thus get closer together for higher energy while the minina get further apart. The two temporal differences are equal at the energy minimum.

All of the interactive graphics presented here provide intuition into a desired parameterization of the exact solution to the pair of nonlinear differential equations.

Uploaded 2018.10.07 — Updated 2018.10.11 analyticphysics.com